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Abstract. We present an iterative method to obtain approximations to Bessel 
i-C . functions of the first kind J p (x) (p > — 1) via the repeated application of an integral 

operator to an initial seed function fo(x). The class of seed functions fo(x) leading to 
' sets of increasingly accurate approximations f n {x) is considerably large and includes 

& . a- n y polynomial. When the operator is applied once to a polynomial of degree s, it 

yields a polynomial of degree s + 2, and so the iteration of this operator generates 
sets of increasingly better polynomial approximations of increasing degree. We focus 
on the set of polynomial approximations generated from the seed function fo(x) = 1. 
^ ' This set of polynomials is not only useful for the computation of J p (x), but also from 

, a physical point of view, as it describes the long-time decay modes of certain fractional 

| diffusion and diffusion- wave problems. 

(N , 

PACS numbers: 02.30. Gp, 02.30.Mv, 02.30.Tb, 05.40.-a 
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1. Introduction 

Bessel functions play a central role in numerous problems in Science [U EJ EJ HJ E]. 
In particular, Bessel functions of the first kind J p (x) appear in the solution of many 
problems involving wave and heat equations, especially in systems with spherical or 
cylindrical symmetry (2j [31 0] . Such functions have been intensively studied for more 
than three hundred years [5] and a large body of results concerning their properties 
is now available [HI [7J [8]. It is well-known that they can be expressed in differential 
and integral forms, in terms of infinite series, as a solution of differential equations or 
recursive relations, etc. Given their ubiquity, their efficient numerical computation and 
its approximation in terms of simpler functions is an issue of great interest [91 [TOj [HI H2] . 
In particular, the approximation of special functions in terms of polynomials has some 
key advantages, as the latter can be effortlessly evaluated and manipulated. 

In what follows we show how to derive infinite sets of polynomial approximations 
to J p valid for p > — 1 [and, of course, also for negative integer values of p by virtue 
of the relation J_ p (x) = (—l) p J p (x) for p integer]. Each set is generated from an 
initial polynomial seed function fo(x) via the iteration of an integral operator and it 
eventually converges towards a properly normalized Bessel function of the first kind, 
J p (x) = 2 p p\Jp(z p x)/(zpx) p , where z p denotes the first zero of J p (x). When applied 
to a polynomial approximation of degree s, the integral operator yields an improved 
polynomial estimate of degree s + 2. The set of polynomial approximations Ba^ we 
primarily focus on stems from the possibly simplest seed function fo(x) = 1. However, 
we also present some results for another set of polynomials Be^ corresponding to the 
choice fo(x) = 1 — x. 

Our method is inspired by the interesting properties of the integral operator 

K[f] = 4l l -^rdw^f{v) (1) 



p J x U 2p+1 J0 

recently introduced by the authors in collaboration with Borrego [13] to study Fourier- 
Bessel solutions of fractional diffusion equations. The operator A p is closely related to 
Bessel's differential equation and it has the remarkable property of leaving the function 
J p (x) invariant (cf. Appendix A). Interestingly enough, when the operator A p is properly 



normalized, the resulting functions converge to J p (x). In the language of dynamical 
systems theory one could therefore respectively speak of "fixed point of the application 
defined by the operator" and "basin of attraction" when referring to J p (x) and the set 
of seed functions which eventually converge to J p (x). 

As an aside, it is interesting to note that our method is similar in spirit to 
Neumann's method for tackling integral equations, where the integral operator defining 
the integral equation is used to generate a function series starting from a zero-th order 
approximation function [To] . In contrast, our method differs significantly from the 
techniques used in recent works on polynomial approximations for Bessel functions. 
For example, Gross [10] exploited the similarity of an integral related to a problem of 
electromagnetic scattering in a conducting strip grating with an integral representation 
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of J p (x) to devise a polynomial approximation valid for Jo(x) and J\{x). Millane and 
Eads subsequently extended this type of approximation to any J p (x) of integer order 
p [TT]. More recently, Li et al. generalized these results by computing approximations 
valid for any real p [12]. At a later stage, we shall use the results in [12] as a reference 
to test the accuracy of our own approximation. 

The present work is organized as follows. Section [2] is devoted to the introduction of 
some preliminary definitions. In section [31 the operator A p is shown to have an attractor 
proportional to J p (x). We demonstrate that repeatedly applying A p to any non-zero 
polynomial seed function generates a set of polynomials which converge to the attractor. 
In section H] we show that the attractor of a slightly modified operator A p is J p (x) itself. 
We subsequently study the set of polynomial approximations generated by two different 
seed functions and discuss some numerical results (section [5]). As shown in section [6l 
one of these sets appears naturally in the context of some fractional diffusion/diffusion- 
wave problems, thereby describing the long-time decay of the solutions. Finally, we 
summarize our main conclusions and briefly outline some avenues for further research 
in section 



2. Preliminary definitions 

For the purpose of finding polynomial approximations to J p {x), it is convenient to 
introduce the following set of normalized functions 

Jp,ni%) 2^p!(zp n x) ^ J p (z Ptn x) , (2) 



where z p>n is the n-th zero of J p (x). The function J p ,i(x) will hereafter play a central 
role, and hence we shall use the short-hand notation J p (x) to denote it. Likewise we 
shall In terms of the above notation, (TSl) can be rewritten as 

x p _, 

J P ( X ) = 7^\ J p( x / z p) (3) 

for n = 1. The last equation implies that, if polynomial estimates for J p (x) are available, 
they immediately lead to approximations of the same type for J p (x). From the series 
representation of J p (x), 

X P oo ( — l) fc fx\ 2k 

W = ^E kl{k + p)l [j) > ( 4 ) 

one easily sees that J p (0) = 1. 

Let us introduce the following linear integral operator 

Ap,»[/]=4„/' J^rf dvv^f(v) (5) 

acting on the set of functions C whose elements are those functions for which the double 
integral performed by A p n is well defined. It is easy to show that any function which is 
equally singular or less singular than x a with a > max[— 2, —2 — 2p] is an element of C. 
[We say that g(x) and x a are equally singular if g(x)/x a —> const, as x — > 0; g(x) is less 
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singular than x a if g(x)/x a — > 0]. Note also that any polynomial function is 

mapped onto another polynomial function by the operator A p n . 



As shown in section Appendix A , the functions J p , n (x) remain invariant under the 



action of the operator A p n , i.e., 



Jp,n{%) Jp,n{%). (6) 



In what follows we shall use the short-hand notation A p = A p i for simplicity. With this 
notation, one has 



A p 



Jp[X j 



J p (x) (7) 
for the special case n — 1. 

3. The operator A p has an attractor proportional to J p (x) 

We aim to show that 

Urn A;[P 8 (x)\ oc J p (x), (8) 

where P s (x) = Yf r =o c r% r is a polynomial of degree s. To accomplish this task, we shall 
first show that Al^l] has an attractor proportional to J p (x). Next, we shall derive an 
expression for A™[P s (x)] in terms of A?[l], . . . , A"[l]. Finally, we shall argue that the 
limit n — > oo of that expression yields (IE])- 

For the sake of simplicity, we shall use the notation A™[1] = I m (x) in what follows 
(in this short-hand notation, no explicit reference to the index p is made, which should 
be inferred from the context). Clearly, the J m (x)'s are polynomial functions of degree 
2m in x. In particular, one of course has A°[l] = Io(x) = 1. In order to study the 
behaviour of these functions as m — > oo, we now invoke a relation obtained from (B4), 
(9) and (11) in reference [13], namely 

\ 2m p— 1 

I m (x) = A-[l] = E ( — ) -TJ—^ -U*) ( 9 ) 




(for a short demonstration, see Appendix B ). Note that, as m gets larger, the weight of 
the large-/c terms in the above sum decreases very rapidly due to the prefactor z~^ m . In 
the limit m — > oo, the only relevant contribution corresponds to the k = 1 term. Hence, 
one has 

I m { x ) (p Jp( x ), as m — * OO) (10) 
where the number ( p is 

1 z p ~ l 

C =— JL (11) 

Because of the linearity of A™ formula fflQ|) can be rewritten as 

rSSo K{ 1 /Cp] = jp(x)- (12) 

Equations (fl2|) and (J7j) provide two remarkable ways of defining p-th order Bessel 
functions of the first kind by means of the integral operator A p : (i) as the limit of the 
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feedback process f n+ i(x) = A p [/ n (x)] with fo(x) = 1/Cp, and (ii) as the solution of the 
fixed point equation f(x) = A p [f(x)]. These definitions nicely resemble the way in which 
fractals are defined in terms of the Hutchinson operator [H] , Of course, equation ( flOl) 
also allows one to compute increasingly accurate estimates for the normalized functions 
J p (x) and thus for J p (x) by virtue of (J3j). 

Evaluation of A%[P s (x)} 
A straightforward application of the operator A" to the monomial x r gives 

AJx r ] = r- r = a r — a r x r+2 for p > —1, r > 0, (13) 

pl J (2 + r)(2 + 2p + r) y ' ~ ' V ' 

where a r = z 2 /\(2 + r)(2 + 2p + r)]. It can be proven by induction that 

A> r ] = EC-l^&fcCr.jOWa;) + (-l) n b n (r,p)x r+2n . (14) 
fc=i 

In the above expression we have introduced 

fc_1 2k r" (2p + r)" 

6 fc (r,p) - n a r+2m = , p (T , + 2fc)! , (r + 2p + 2fc)!! (15) 

as well as the definitions (2n + 1)!! = (2n + l)!/(2 n n!) and (2n)!! = 2 n n!. Note that 
b n (r,p) goes rapidly to zero for large n: 

yr,p)~( 2p /2) 2 7N) 2 (16) 

Applying A™ to an s degree polynomial P s (x) = J2r=o c r xT an d using (fl5l) and ( IT41 . one 

gets 

A£[P s (x)] = E(-l) fc - 1 /„- fc (x)Ecr6fc(^P) + (-l)V n E c A(r,p)x r . (17) 

fc=l r=0 r=0 

Next, let us investigate the behaviour of the above expression when n — > oo. First 
note that, according to ( JIB]) , the coefficients b n (r,p) roughly go to zero as {n\)~ 2 for 
large n. As a result of this, the second term on the r.h.s. of ( TT71) becomes negligible 
with respect to the first one, which is a polynomial of degree 2n — 2 with an independent 
term (of course, the larger the value of x, the larger the value of n necessary to make 
this term negligible). On the other hand, it is not difficult to see that the first term 
of (jT7l) tends to an expression proportional to J p (x), since (i) according to ( JTUI) . when 
n — k — > oo the polynomials I n -k(x) tend to the limiting expression ( P J P , and (ii) those 
terms in the sum for which n — k is small, i.e., those terms for which k is large, become 
negligible due to the fast decay of bk(r,p) with increasing k. We thus conclude that 
Ap[P a (a;)] cx J p (x) when n — > oo. 
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4. Generation of J p (x) via the normalized operator A p 

As we have just seen, when successively applied to a polynomial, the operator A p 
generates another polynomial which rapidly approaches J p (x) up to a prefactor. A 
possible way to get rid of the prefactor is to divide these polynomials by their value 
at the origin. The resulting "normalized" polynomials becomes equal to one at x = 0, 
which is precisely the value taken by J p (x) at x = 0. This prompts us to introduce the 
"normalized" operator 

kM = (18) 

acting on the subset of functions f £ C for which A p [f] x=0 (that is, A p [/] evaluated at 
x = 0) is non-zero. In particular, using ( ITfl) . we have 

in [p M1 _ ELi(-l) fc " 1 ^- fc (x)S^ c r 6 fc (r,p) + (-l)"x 2 "E; =0 cA(r, P K 

and Ap[P s (x)] — > J p (x) as n — > oo. In view of the above derivation, one expects that 
any well-behaved f(x) (in the sense that it can be approximated arbitrarily well by a 
polynomial) also converges to the same attractor, i.e., 

A£[/(x)] -> J P (x), n^oo. (20) 

The iteration of the operator A p will allow us to generate families of approximations 
to J p (x) by using different seed functions fo(x). Each of the seed functions fo(x) leads 
to a series of functions {fo, fi, f 2 , ■ ■ .}, with /„ = A"[/ ], that converges to J p (x) as 
n — > oo. One expects that the convergence to the attractor becomes faster as the initial 
condition gets closer to J p (x). On the other hand, from a practical point of view it would 
be desirable that the chosen seed function / is simple enough to ensure that the integrals 
resulting from the iteration of A p are known and easy to calculate analytically. Obvious 
candidates are polynomials. Surely, the most simple initial function is fo(x) = 1. This 
function gives rise to the set {f n (x) = Ba^(x)}, defined as 

Ba^(x)^A p «)[l] = ^||. (21) 

The initial condition fo(x) = P s (x) = 1 is a particular case of a polynomial function for 
which c = 1 and c r = (r > 0). Inserting this expression into (jT7j) and recalling the 
definition of I n (x) one gets 

n 

In(x) = ^-\) k - x b k ^,p)I n _ k {x) + (-l) n b n (0,p)x 2n , (22) 

k=l 

where bk(0,p) = Zp k p\/[2 2h k\(p + k)\] [cf. fTl5|) ]. Inserting this expression into ( 122"]) one 
finds 

= U~ iri ^£^- In - t(x) + { - ir ^£^r x2n - (23) 
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The first few polynomials Ba^(x) computed from (1211) and (1231) read 



o \ 

(p)(„\ _ i _ ji 



Ba^ ; (x) = 1 



Bar(^) = 1 - x 
Ba« W = 1 - + £±1*«, 



As a test of the present method, we shall also investigate the behaviour of another 
set of functions Be^(x) generated from the fo(x) = 1 — x. In this case, from (|T7|) and 



the definition of A™ one finds 



A n n , = SLi(-l) fc "^n- fc (x)[6 fc (0,p) - 6 fc (l,p)] + (-l) n x 2n [b k (0,p) - 6 fc (l,p)s] 

Pl ' J__ U=l(-l) k - 1 In-k(0)[b k (0,p)-b h (l,p)\ { ] 



The corresponding polynomials can be easily evaluated via the equations (1231) and (I15j) . 
The first few polynomials are given below: 

Bej^O) = 1 -x, 

Bei^) = l-^* 2 + fc^, 

Be(p) = _ 10(p + 2)(2p + 5) 2 2 5(p + l) (4p 2 + 16p+15) ^ 
2lJ 3(V + 36p 2 + 115j9+113) 4p 3 + 36p 2 + 115p + 113 

32(p+l) 2 (p + 2) 5 

3v • 



3(4p 3 + 36p 2 + 115p+113) 
5. Numerical results 

Li, Li and Gross (LLG) proposed in Ref. [12] an approximation to J p (x) based on its 
integral representation. They obtained an infinite series whose truncation leads to the 
following polynomial approximation of J p (x), 

T(P) = Vf If n^m + n-l)! /x\ 2 ^ 
n ^ ] m\(n-m)\T(m + p+l) \2j { ] 

thereby providing an alternative approximation to truncated Taylor series. Figure 1 in 
[12] depicts a comparison of J p (x) with the LLG polynomial approximation of degree 
2 x 10 +p and with the Taylor series truncated to the same degree for p = 0, 3/2, 3, 5. In 
figure [T]of the present paper we have superposed their results with our own results based 
on the approximation obtained by making the replacement J p (x/z p ) — > Baio }(x/z p ) in 
Eq. ([3]). We see that our approximation performs very well over a larger range of x 
values than the LLG polynomial approximation and the truncated Taylor series of the 
same order. Even though the LLG polynomial approximation oscillates around J p (x) 
over a larger x interval than the other two approximations, it starts deviating from the 
exact curve at smaller values of x. In all cases, deviations from J p (x) are shifted to larger 
x values with increasing order p. In this sense, all three approximations become better 
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X 

Figure 1. Comparison of J p (x) (solid line) with the LLG polynomial approximation of 
order 2 x 10 +p [12], L±q, (short dashed line), the truncated Taylor series (dotted line) 
and the polynomial approximation Ba^ 10 ^(a;) of J p (x) (dashed line) [see (|2Tj) and (|23l) ] 
for p — 0, 3/2, 3, 5. In the lower panel the line corresponding to our approximation for 
J§{x) lies on top of the exact one. 

with increasing p. Regarding the computational efficiency, the CPU time employed for 
calculating the polynomial approximations obviously depends on their order n. For 
example, when using the program Mathematical the CPU time necessary for computing 
the curves in the top panel of figure [1] (corresponding to n — 10) is roughly the same 
for each of the three polynomial approximations depicted; this time is roughly one half 
of the time required to compute the curves corresponding to n = 20 and, also, one half 
of the time required when using directly the Bessel function algorithm. 

In figure |2] [figure [3] we compare Jq{x) [^(z)] with the polynomials Ba^x) 
[Ba^-^x)] for n = 1,2,3,4,5,10. For fixed values of n and x the approximation is 
more accurate for p = than for p = 2. For other p values we have checked that, 
in general, the accuracy of our approximation for J p (x) decreases with increasing p, 
as opposed to the behaviour observed for J p (x) (see figure [p. Clearly, the different 
behaviour is related to the rescaling introduced by (|3]). 

Finally, we proceed to compare J\ (x) with the polynomial approximations Ba^ (x) 
and Be^\x) for n = 1,2,3,4,5,10. This comparison is shown in figure HI Not 
surprisingly, the family Be^\x) performs better than Ba^^x) since it starts from the 
seed function f (x) = l — x which is a better approximation to Ji(x) than fo(x) = 1. 
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Figure 2. Comparison between Jo{x) (solid line) and Ba^(i) with n = 1, 2, 3, 4, 5, 10. 
The value of n corresponding to each line is shown. 




Figure 3. Comparison between ^(a;) (solid line) and Ba^(i) with n = 1, 2, 3, 4, 5, 10. 
The value of n corresponding to each line is shown. 




Figure 4. Comparison between J±(x) (thick solid line) and Ha,^\x) (dashed lines) 
and Be^ -* (x) (thin solid lines) for n — 1, 2, 3, 4, 5. The value of n corresponding to each 
line is shown. 



6. The polynomials Ba^(x) describe fractional diffusive and oscillatory 
modes 

As shown in [13] . the polynomials Ba„ (x) can be used to express the long-time decay 
modes in some fractional diffusion problems in a d- dimensional sphere. In what follows 
we shall demonstrate that the Ba^(x)'s also appear in problems whose solution is given 
by the fractional diffusion-wave equation subject to the same geometry and boundary 
conditions (note, however, that for diffusion-wave problems one must additionally specify 
the initial velocity). Take, for instance, the problem described by the equation [16] 

^^ = ^ 7 V 2 c(r,t), (26) 

where K 1 is a coefficient and the operator (P/dt 1 denotes the Caputo fractional 
derivative of order 7 [161 El HH] : 

— ; — ^7 c / 7 r — — ; — dr. n — 1 < 7 < n (27) 

dP T(n - 7) Jo (t - T p +1 ~ n ' ' y J 

with n an integer. In [13], problem B in Section III, the long-time solution of (|26[) 
corresponding to the boundary condition c(R, t) = and the initial condition c(r, 0) = c 
was discussed for the case < 7 < 1. Physically, c(r,t) represents the decay of 
a homogeneous initial particle concentration inside a hyperspherical volume with an 
absorbing boundary of radius R. The solution consists of a series of decay modes whose 
spatial part is expressible in terms of the polynomials Ba^(x) via a recursion relation 
essentially identical with (1231) . 
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In the range of values 1 < 7 < 2, (1251) becomes a diffusion-wave equation 
P31 H21 123 I2U- For c(r, 0) = c and dc(r, t)/dt = at t = the solution reads 

c(r, t)/c = 2(r/i2)"' £ ^ ^ 2S 7 [-(^/i?) 2 *^] (28) 

j=l Z rj,j J rj+l\ z r},j) 

with 77 = d/2 — 1, which is also the solution for < 7 < 1 (fractional diffusion equation) 
when c(r, 0) = . The solution is obtained by separation of variables and subsequent 
use of the relation [16], [T7] 

— £ 7 [-u/^] = -uFE, [-u 7 fr\ (29) 

for the derivative of the Mittag-LefHer function E^[-\. Inserting the asymptotic expansion 
of the Mittag Leffler functions 

00 / 1 \m+l 

^H-E r l (1 _ mT) ^ (so) 

for < 7 < 1 and 1 < 7 < 2 into ( 1281) and grouping the terms with the same power of 
t, one finds 

eMV^fP^^Vw*) (3D 

for < 7 < 1 and 1 < 7 < 2 (in the above equation I m (0) must be evaluated for 
p = if). Thus, the spatial dependence of the long-time solution can be expressed in terms 
of a suitable superposition of the "fractional modes" Ba^(x). The above fractional 
solution extends the classical solution (cases 7 = 1,2), where the spatial modes are 
proportional to Bessel functions. In other words, the role of the m-th d- dimensional 
normal mode J v (z v ^ m r / R) / (r / R) v in the normal diffusion equation (7 = 1) and wave 
equation (7 = 2) is played by the m-th rf-dimensional fractional mode Ba^ (r/R) in 
the fractional diffusion equation (0 < 7 < 1) and fractional diffusion- wave equation 
(1< 7 <2) : 

Jvi ^ v R ' ] ^B^(r/R). (32) 
7. Summary and Outlook 

The present work deals with a novel method to obtain polynomial approximations to 
Bessel functions of the first kind. In some cases the obtained polynomials turn out to 
be more accurate than truncated Taylor series and the polynomials in [12] over a wide 
range of parameter values. We have seen that the polynomials Ba^ (x) generated by the 
initial function fo(x) = 1 are interesting in their own right, as they represent the spatial 
modes describing the long-time behaviour of certain solutions of fractional diffusion and 
diffusion-wave equations. 

From a computational point of view, a nice property of our method is that for fixed 
x and p the distance Ba^ (x) — Bajf •* (x) | between successive approximations decreases 
rapidly with increasing iteration number n, and one can set a threshold value below 
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which, for practical purposes, convergence to the corresponding Bessel function may be 
considered to have taken place. 

Integral operator methods developed along similar lines might be able to provide 
alternative polynomial approximations in the range p < — 1. Of course our solution is 
also valid for negative integer values of p by virtue of the relation J- P (x) = (—l) p J p (x) 
for p integer. Such techniques might also be useful for other kinds of Bessel functions 
(e.g. Bessel functions of the second and the third kind). An open question is whether 
similar iterative methods can be applied to generate polynomial approximations to other 
special functions associated with ordinary differential equations. 
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Appendix A. Invariance of J p n (x) under the action of A p n 

The starting point is the differential equation fulfilled by the Bessel functions 



y 



dy 2 dy 

Using the transformation y = z n , p v in the above equation gives 



i d J<n( Zn_n. 



+ v- 



djv%\ Zn.n V ) 



dv 2 dv 
Multiplying by v p ~ x jz 2 n and rearranging terms we get 



+ (4,n v2 -P 2 ) J P ( Z P,nV) = 0. 



p+1 d 2 Jp(z p , n v) v p dJ p (z p , n v) p 2 v p 1 



dv 2 



+ 



p,n P,n 

The above equation can be rewritten as follows 



dv 



Jp(.Zp,n V) ■ 



p.n 



V Jp^Z pn V ) 



1 d 



Z p,n dv 



V ^ Jp{,Zp )U v) 



dv 



Integrating between and u we get 



2T 

U 2p+l Ji> 



d r 



dw P+1 J p (z Pin v) = — 



where we have used 



\imv 2p+l ^- (v- p JJz pn v)) = lim 
^->0 dv ^ ' v-yO 



U Jp{^Zp n tt) 

1 dJ p (z n v) 



dv 



- pv p J p (z Ptn v) 



z P)n lim v p+l J p+1 (z Pin v) = (p > -1) . 



Integrating once again (IA.5j) between 1 and x and using J p (z Ptn ) = we get 
- 1 du 



lx U 2p+1 JO 

Finally, multiplying the above equation by 2 P p\z~ p yields ((6]). 



(A.l) 



(A.2) 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



Appendix B. Proof of formula ([9]) 

The Fourier-Bessel expansion of a function g(x) is given by g(x) 

2 n 



Ck 



Jp+l( z p,k) JO 



xg(x)Jp(z p ^x)dx. 



Efcli c k Jp{zp, k x) with 
(B.l) 
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Therefore, OH]) is equivalent to the statement that the coefficients of the Fourier-Bessel 
expansion of x p I n (x) are 

/ \2n p-1 

c k = c(k,p tn ) = 2[-^) _ (B.2) 

\ z P,k / Jp+l\ z p,k) 

We are going to prove this equation by induction. To start with, it is well-known [6] 
that the coefficients of the Fourier-Bessel expansion of x p Iq(x) = x p are 

c(fc,p,0)= (B.3) 

This justifies ( IB. 21) for n = 0. Therefore, to prove equation (jHJ) is equivalent to prove 
that c(k,p,n + 1) = (z p /z Pt k) 2 c(k,p,n). By using (IB. II) with p(x) = x p I n+ i(x) and 
making the substitution y = z p ^x inside the integral one sees that 

2 f z p,k 

c(k,p,n + l)= p+2 / y p+1 J p (y)I n+1 (y/z Ptk )dy. (B.4) 

Z p,k J p+A z P,k) J ° 

Using y p+1 J p (y) = d[y p+1 J p+ i(y)]/dy, and integrating by parts one gets 

2 r z v,k d 

c(k,p,n + l) = — 2 -/ y p+1 J p+1 (y)—I n+1 (y/z p , k )dy (B.5) 

as the boundary terms vanish. But, from (jSJ), 

W (x) = 4l l ^Ti I" dvv 2 *+%(v). (B.6) 



u 2p +K .. 

Inserting this expression into (IB. 51) . using the relation J p+1 (y)/y p = —d[J p (y)/y p ]/dy 
and integrating by parts one gets 



2 



c(k,p,n + 1) = ^ fc / y p+1 J p (y)I n (y/z Pik )dy (B.7) 

J p+A z P,k) Jo 

as the boundary terms vanish. Comparing this result with the expression of c(k,p,n) 
given by (IB .41) one sees that c(k,p,n + 1) = (zl/z pk )c(k,p,n), which is just the result 
we aimed to prove. 



